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Abstract 

We study how bubbles grow after the initial nucleation event in generic first-order 
cosmological phase transitions characterised by the values of latent heat L, inter- 
face tension a and correlation length £, and driven by a scalar order parameter (p. 
Equations coupling 0(t, x) and the fluid variables v(t, x), T(t,x) and depending 
on a dissipative constant T are derived and solved numerically in the 1+1 dimen- 
sional case starting from a slightly deformed critical bubble configuration 0(0, x). 
Parameters L, a, £ corresponding to QCD and electroweak phase transitions are 
chosen and the whole history of the bubble with formation of combustion and 
shock fronts is computed as a function of T. Both deflagrations and detonations 
can appear depending on the values of the parameters. Reheating due to colli- 
sions of bubbles is also computed. 
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1 Introduction 



The bubble nucleation in cosmological first-order phase transitions as well as the prop- 
agation and stability of planar interfaces have been discussed extensively in the litera- 
ture [pH — [1T9| - The purpose of this paper is to join these two ends of the life of a bubble 



by giving a dynamical model which describes the entire history of the bubble from the 
initial configuration via initial acceleration into a large bubble growing with constant 
velocity. 

The stage for the events in this paper is the cosmic fluid with a first-order phase 
transition at T = T c — in practice either the QCD or the electroweak (EW) phase 
transition. For physical quark masses there is no symmetry associated with the former 



one and its order is not definitely confirmed [ 20[ . The latter is a symmetry breaking 



transition and, at least for not too large Higgs masses, it is of first order both on the 
basis of perturbative [PI]1 -|PB| and nonperturbative lattice Monte Carlo PT|-P5| work. 
Inflationary transitions are essentially vacuum ones and the considerations here do not 
apply. 

The main quantity characterizing the cosmic fluid is its energy-momentum tensor 
T^ u = (t+p)u^u v '~pg^ v . The first-order nature implies that there are two phases, a high 
temperature phase (symmetric, quark-gluon plasma phase) with pressure p q (T) and a 
low temperature phase (broken symmetry, hadron phase) with pressure ph(T), which 
can coexist at T c : p q (T c ) = Ph(T c ) but p' q {T c ) > p' h {T c ). We shall describe the transition 
by a scalar order parameter field (j)(t,x). This is obvious for a symmetry breaking 
transition, but we shall use the same description also if no symmetry is involved: the 
order parameter could then be, for example, the energy or entropy density. The bubbles 
are configurations of 0(t,x). 

The problem now is to derive equations of motion for the total cosmic fluid - 
order parameter field system. We carry this out in analogy with the reheating problem 



in inflation |KJ, |3T|. The total energy momentum tensor of the system is conserved 
(d v T^ u = 0) but those of the fluid and order parameter field subsystems are not. 
Physically, entropy produced at the bubble interface couples the behavior of (f) with the 
fluid. The strength of this coupling is described by a dissipative constant V. It has been 
related by a fluctuation-dissipation formula to equilibrium averages in refs. f32|, |33|j, but 
we use T as a phenomenological parameter. Estimates for it in the EW theory have 
been given in refs. |TT| , |34| - 

Given the equations of motion one can solve them numerically and study how 
bubbles corresponding to given initial supercooling (which follows from nucleation anal- 
ysis), given parameters of the transition (latent heat, interface tension, correlation 
lengths), and given dissipative constant V evolve. Collisions of bubbles can be simi- 
larly studied. In this first discussion we shall solve the equations in 1 + 1 (one time, 
one space) dimensions, which contains the main qualitative features. For complete 
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quantitative results one must go to 1+2 or 1+3 dimensions. 

The initial stages of bubble growth have also been computed by integrating 
spherically symmetric 1+3 dimensional hydrodynamic equations in ref. ||35|| . The ve- 
locity is taken as a free parameter, effectively parametrised by the magnitude of energy 
flux. Our work differs in one essential aspect from this: the dissipative constant T is 
at least in principle calculable from the theory. 

The paper is organized as follows. In Section |2| we introduce our model for 
the first-order phase transitions. In Section ^| we review the general hydrodynamic 
conditions for the bubbles, and describe the two kinds of solutions, detonations and 
deflagrations. In Section £| we describe our results for different time-dependent phe- 
nomena like the initial stages of bubble formation, the sharpening of shock fronts and 
the collisions between expanding bubbles. In Section [| we give an account of the differ- 
ent steady-state variables for deflagrations (temperatures and velocities) as a function 
of T. The conclusions are in Section ^. 

2 Equations of motion for the cosmic fluid and the 
order parameter field 

The system we consider contains the cosmic fluid which has supercooled in the metastable 
high T phase (call it by convention q) to some temperature Tf. At this temperature 
nucleation of bubbles of the low T phase (call it by convention h) becomes sufficiently 
frequent for the phase transition to effectively take place. We shall first define the 
quantities appearing in the equations of motion. 

The bubbles are defined as configurations 0(t, x) of a scalar order parameter <p. 
The (meta) stable states of the system are defined by the minima of the effective poten- 
tial V(4>, T) of the order parameter <fi. The equations will be formulated for a general V, 



but for numerical calculations we shall use a quartic parametrisation p, 14 



V(4>, T) = i 7 (T 2 - T o 2 )0 2 - i«T0 3 + iA0 4 . (1) 



The full functional form of the ring summation improved effective potential ||22||— 126} 



deviates somewhat from this, and even more does the effective potential containing 
nonperturbative effects [E]J. The methods developed here can be straightforwardly 
extended to these improved potentials. 

The physical quantities corresponding to the parameters a, 7, A and T in eq. ([]]) 
are the latent heat L, the interface tension a, the critical temperature T c and the 
correlation length £(T C ) = £ c |L5[j . The quartic parametrisation in eq. (p]) implies that 



the correlation lengths in the two phases at T c are equal, but actually they are different 
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The primary physical quantities characterizing the transition are T c , L, a and 
£ C) and given these, we can always solve for the parameters a, 7, A and T in eq. (p. 
In this way we can use eq. ([!]) also for QCD. Its use for a symmetry breaking transition 
is, of course, obvious. 

For the equation of state of the system we shall take 

Pq (T) = aT* Ph (T) = aT* + B(T), (2) 

where B{T) = — V((p m i n ,T) is the difference between the free energy densities of the 
symmetric and the broken symmetry phases. The number of effective degrees of free- 
dom of the high temperature phase is denoted by g*, and a = (it 2 /90)g*. 

Summarizing, our dynamical variables are the four- velocity of the fluid x), 
the scalar field (f>(t, x) and the local temperature T(t, x). We use the notations w r = 
4aT 4 and p r = aT A for the pure radiative enthalpy and pressure, respectively. 

To motivate and to write down the equations of motion we first write the effective 
potential of eq. (P in two parts: 

V(4>, T) = ~7T O 2 2 + \\^ + l -!T^ - i«T0 3 . (3) 

" v ' V v ' 

V (<j>) V X {4>,T) 

The part Vo(</>) is related strictly to the scalar field, whereas the part Vi(<f>,T) includes 
the interaction of the field with the thermal bath. The total energy-momentum tensor 
is 

= d^d^-g^d^d^-VM] (4) 

+ [w r - T dVl ^ T T) ]u^ - grfa - Vi((f), T)] , (5) 
and it is of course conserved: 

d^ v = . (6) 

The idea now is to bring in a dissipative term which acts as a friction force for 
the scalar field and accounts for the entropy production at the phase transition surface. 
This is done by splitting eq. (|6]) into two: 

= 8^ = [d^% + [^THrad = 8»-8 v . (7) 

The form of the Lorentz-covariant dissipative term 5 V is adopted from the context of 



inflation |30|, |3T|]. Because of the temperature dependence of the effective potential, 
the choice of the terms [d^T^]^ and [9 M T M ^] ra d in eq. (0) is not unique. It seems most 
natural to make the splitting in the following way: 

dp{d»<j>dr4> - g ^[\d a <pd a <p - v Q ((f>)]} + = -^d^d^ 
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After all, if Vi(4>, T) could have been written in the form Vi((f>,T) = V^(0) + Vr(T), 
there would be no ambiguity in the splitting, and exactly eqs. (||) would result. Our 
final equations are then 

W + f = -^d,<p 

d,{[w r - T%]u»u» - g^[p r - V}} = (±it"3^ + %)d v (j> , 

where a common factor has been dropped. 

Contracting both sides of the lower of eqs. (j|) with the fluid four- velocity u v one 

gets 

T ^[(^-|^K] = ^K^) 2 , (io) 

where s r = 4aT 3 is the radiative entropy and — |^ the entropy associated with the 
order parameter. This equation relates entropy production and the gradients of <p via 
the constant T. Note that in a weak coupling theory ]TT], Q 1/T ~ l/r c ~ nva ~ g 2 T. 

In this paper we study eqs. @ in 1+1 dimensions, which corresponds to pla- 
nar symmetry. While this is a drastic simplification, it nevertheless should correctly 
describe the late stages of the bubble growth in the 1+3 dimensional world. Planar 
symmetry also allows us to compare our results with analytical calculations. 

In the planar-symmetry case it is also illuminating to write down the equations 
for the steady-state solution. At large times the system should evolve to a solution 
containing a combustion front moving at constant velocity. In the rest frame of the 
front all time derivatives then vanish and eqs. (|9|) become 

4f(*) = ^ + (ii) 

dV 

(4aT 4 -T—) 1 2 v = const. (12) 

(4aT 4 - T^h 2 v 2 + aT 4 + ^'(x) 2 - V = const. (13) 

Eqs. ( |i~2"D and (0) actually also follow from the steady-state energy-momentum con- 
servation: 

d x T x ^ = . (14) 
Similarly, the entropy production equation ([[(]) becomes 

T±[(s r - ^)jv) = ^ 2 v 2 W(x)] 2 . (15) 

The standard analysis of deflagration and detonation bubbles is a study of 
what solutions eqs. (0), flHf ) allow. For given initial Tf this leaves a one parameter 
family of solutions. For detonations there are four quantities, T q , v g , T h , Vh, constrained 
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by two equations and by the boundary value T q = Tf. For deflagrations the shock 
front has to be taken into account, and there are two more quantities but also two 
more equations. This is discussed in Section |3|. Eq. (|ll|) is the new one which gives the 
new physics permitting one to choose the correct solution within the one parameter 
family. 

3 Deflagrations vs. detonations 

There are two kinds of bubbles allowed by the hydrodynamics. These are called de- 
flagration and detonation bubbles |?[] according to the nature of the phase transition 
front. 

Consider fluid flow in the rest frame of the phase transition front. The incoming 
flow velocity is denoted by v\ and the outgoing velocity is denoted by v 2 (see fig. 1). 
In a deflagration the incoming flow is subsonic, v\ < c s , and the fluid is accelerated by 
the phase transition, u 2 > V\, whereas for a detonation the opposite is true: v\ > c s 
and t>2 <V\. Depending on whether the outflow is sub- or supersonic, these processes 
are further divided into weak (t> 2 < c s ), Jouguet (t> 2 = c a ), and strong (t> 2 > c s ) defla- 
grations, and strong (t> 2 < c s ), Jouguet (t> 2 = c s ), and weak (t> 2 > c s ) detonations 0. 

Consider then the structure of the bubble in the rest frame of the ambient fluid. 
When the bubble has grown large enough, any memory of the initial shape of the 
nucleated bubble should be lost. The bubble can then be described as a similarity 
solution of the hydrodynamical equations, i.e., it expands linearly with time, otherwise 
maintaining its shape and profile. The fluid has to be at rest both at the center of 
the bubble, and far away. Thus, in a deflagration bubble, the phase transition front is 
preceded by a shock wave which heats up the fluid and sets it moving outward. The 
phase transition front then brings the fluid back to rest (see fig. 2). In a detonation 
bubble, the fluid is at rest when it is hit by the phase transition front, which leaves 
the fluid flowing outward. A rarefaction wave follows, bringing the fluid at rest (see 
fig. 3). We denote the velocity of the phase transition front in this frame by t>d e f for 
deflagrations and by Vdet for detonations. Weak deflagrations have t>d e f < c s , whereas 
for strong deflagrations v^et > c s . Detonations always have t>d e t > c s . We can exclude 
strong detonations ||, because they leave the fluid flowing too fast: no such similarity 
flow exists that would bring the fluid at rest at the center of the bubble |§. For a 
deflagration bubble, the velocity t> s hock of the shock front is also of interest. In 1+1 
dimensions the fluid flows at a constant velocity t>fl U id between the shock and phase 
transition fronts. 

For deflagrations, the velocities t>def, ^fluid and f shock are related to the velocities 
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v\ and v 2 by the equations 



^shock ^'l | shock ^def ^2 1 phase wall 

V1—V2 



ffluid 



l-ui«a shock Vflmd l-v 2 vi phase wall 



where the words "shock" and "phase wall" indicate the front at which v± and v 2 are 
measured. For the phase transition front, we sometimes use the notations V\ = v q 
and t> 2 = Vh, and T q and T h for the temperatures of the incoming and outgoing fluids, 
respectively. The temperature between the fronts is then T q and the temperature of 
the hadron phase is Th (see fig. 2). 

The initial condition, matter at rest in the q phase, at temperature Tf < T c , 
and the equations of state of both phases, do not fix the rate of bubble growth (v dc{ or 
Viet), or the temperature inside the bubble. For simplicity, we illustrate this with the 
bag equation of stated 

Ph (T) = a h T* + L/4 p q {T) = a q T i 

(16) 

e h (T) = 3a h T 4 - L/4 e q (T) = 3 S T 4 
which one gets from eq. (0) by making the small-supercooling approximation 

BiT) ~ t fl - (17) 



Then a q = a, = a — L/4T C 4 . 

The conservation of energy and momentum, and the non-negative entropy pro- 
duction at the phase transition front restrict the possible values of the incoming (ei) 
and outgoing (€2) energy densities (see fig. 4). Detonations require a certain amount 
of supercooling. If the latent heat L is large, the required supercooling can be quite 
substantial, and the h matter then at a highly superheated state immediately behind 
the phase transition front. This has led to the conclusion that deflagrations are the 
more likely process in the QCD phase transition in the early universe [0. 

However, if the latent heat is small, detonations require less supercooling. The 



nucleation temperature can be estimated [15| from 



T t A 



T c J 171 -41n(171 3 / 2 /A) 



where 171 = 41n(t c T c ) and A = -y/l67r/3(cr 3 / 2 /Lv / 7^). The values of the parameters 
a and L are essentially unknown, although results from lattice calculations can be 



lT Jsually the bag constant B — L/4 appears on the q side, with the opposite sign. This normal- 
ization of the zero-point of energy does not affect the hydrodynamics. 
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employed in giving rough estimates for them. In fig. 5 we show the region in the (L, a) 
parameter plane, where weak or Jouguet detonations would be allowed. 

To choose among the allowed solutions (see fig. 4) the internal mechanism of the 
phase transition front needs to be considered. This is the purpose of our model pre- 
sented in Section |2], with the additional parameter T, which will pick a single solution. 
In the next sections we turn to numerical results obtained for this model. 

For the QCD phase transition we used the "qh" parameters 



L = 2 o = 0.1 
f c = 1 5* = 51.25 

Here and in the following all quantities are expressed scaled with powers of T c to make 
them dimensionless: L = L/T^,a = cr/T c 3 ,£ c = £ C T C . From eq. (pT8[), these parameters 
correspond to the nucleation temperature Tf = 0.9943. The values of the latent heat 
L and the surface tension a are suggested by pure glue lattice Monte Carlo simula- 
tions |S^, |37], |38|. The parameter £ c shows up neither when the nucleation temperature 



is calculated (in the thin wall approximation, eq. fllq) ) nor when the steady-state vari- 



ables are calculated (in eqs. (|30H35|)). However, it determines the thickness of the 
phase transition surface. Because the transition is only weakly first order, the actual 
£ c might be larger than our value. Notice that in addition to strongly interacting de- 
grees of freedom, the parameter includes weakly and electromagnetically interacting 
degrees of freedom. Because the mean free paths of these particles are much larger 
than those of strongly interacting particles, these degrees of freedom are actually not 
active during the early stages of the phase transition [^, |39|. However, since we are 
mostly interested in the final stationary stages of the phase transition, all the degrees 
of freedom are included. For these parameters we expect deflagrations only (see fig. 5). 

In the EW case, we assumed that aw ~ 1/30, m top = rnw and tuh ~ 40 GeV. 
The small effective Higgs mass is necessary in order to allow for a generation of the 
baryon asymmetry. Using the improved effective potential we get our "ew" pa- 
rameters 

a = 0.0162 7 = 0.1309 
A = 0.0131 0, = 106.75 



This corresponds to the nucleation temperature Tf = 0.9957 JT2|. 

For the above qh parameters the lowest temperature To where the symmetric 
minimum = still exists is T = 0.8771. In the ew case we have T = 0.9828. That 
these numbers are not too low indicates that our use of the quartic effective potential 
V{<f>,T) is justified [0 eq. (2.21)]. 
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4 Time-dependent phenomena 



To integrate eqs. (§) we wrote a simple 1+1 dimensional relativistic hydrodynamics 
code following Wilson [EJ . Thus we use explicit differencing with operator splitting 



for the hydrodynamic equations. The code variables are 0, 7r = dt<p, E = 7[3aT 4 + 
V((j),T) — T|^t], and Z = 7 2 t>[4aT 4 — T|^]. The velocity t> is solved from E and Z. 
The temperature T (for each grid point) is solved from E and (f) using the functional 
form of V(4>, T). This value for T is then used in V(4>, T) for evolving <fi. The transport 
terms for E and Z are handled with the FCT method f£|. 

We use reflective boundary conditions (see fig. 6). The center of the initial 
"bubble" of new phase is placed at one end of the grid. Allowing the moving wall 
to reach the other end simulates the collision with another similar bubble. The code 
corresponds to a planar geometry, so these are not true spherical bubbles. In this paper 
we study the motion of a planar phase wall. 

4.1 Initial conditions 

Before starting the actual integration, the initial conditions have to be specified. The 
initial bubble has to be larger than critical to start growing, but it is not clear to what 
extent there is a fluctuation in the temperature associated with the fluctuation in the 
order parameter. Possibly the temperature is a bit higher near the critical bubble than 
farther away from it, because latent heat is released in the formation of the critical 
bubble. We can estimate typical temperature fluctuations with classical fluctuation 
theory fj3| : for the quark matter equation of state 



W> = _JL_ (19) 

T 2 YlaTW ' 1 } 

For the qh case, with a radius R c ~ 5£ c of the critical bubble, we get AT « 0.005 
in units of T c . Comparing with Tf in the previous section, this is very large. For the 
ew case we get AT rs 0.001 [p8||. However, we think that it is most straightforward 



and in the spirit of the nucleation calculation (e.g. in ref. [12]) to assume that the 
temperature inside the initial bubble is just Tf. What is most important, the details 
of the initial bubble have no effect on the final steady-state configuration and the 
asymptotic variables T q , Th, f shock, vamd and t>d e f ( or ^det), if only the nucleated bubble 
starts growing. 

Assuming as initial conditions that the fluid velocity vanishes everywhere and 
that the temperature is constant and equal to Tf, we still have to decide the shape 
and size of the initial bubble. These variables have some significance during the early 
stages of bubble evolution, since they affect the initial shape of the shock front and 
also its initial acceleration. In 1 + 1 dimensions it is possible to analytically find the 
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extremum bubble of the effective action by solving the equation (j>"(x) = d^V^^Tf). 
With M 2 = 7 (T| - T 2 ), 5 = aT f and A = 9XM 2 /25 2 the solution is 



where f(x) = (1 — Acoth 2 Mx)/(l — coth 2 Mx). A bubble obtained from this by 
increasing both the amplitude and radius by a small factor (by 5%) was normally 
used as the initial configuration. We experimented also with the exact extremum 
bubble and with a bubble smaller than this one. The exact extremum bubble did not 
evolve anywhere and the subcritical bubble collapsed, leaving behind a disturbance 
in the temperature and flow velocity propagating outward with sound velocity. This 
disturbance is caused by the fact that matter has to flow inward to fill in the area of 
lower energy density from where the subcritical bubble has disappeared. 

4.2 Numerical results 

In fig. 7 the initial stages of bubble formation are shown for the qh parameters. Imme- 
diately after the nucleation, a shock front is originated which spreads out information 
about the nucleation. At first the shock front is not sharp. The temperature starts to 
rise inside the bubble and very soon the phase transition front begins to get shape. At 
about the moment t = 100(1/T C ), both the phase transition front and the shock front 
are clearly visible. 

To get a more precise picture of the phase transition surface, we can use eqs. (|TT|— 
P~3|) . (The solution of these eqs. is discussed in Section @). In fig. 8 the order parameter 
<f)(x), the temperature T(x), the velocity v(x) and the quantity drV are shown in the 
rest frame of the phase transition surface. Hadron phase is on the left and quark phase 
on the right. The width of the surface layer is a few correlation lengths. The curves 
resemble the tanh(:r) -function but for some other parameters the resemblance is not 
as clear. Specifically, the temperature and velocity distributions lose their symmetry 
and are shallower on the hadron side where the effective potential V(<p, T) is non-zero. 
The "center points" of these distributions are not quite at the same place as that of 
the order parameter, but are shifted towards the hadron phase. For some parameters 
the shift can be of the order of £ c . 

In fig. 9 the development of the shock front is illustrated. The shock front 
is shown at times t = 160 and t = 3840. At early times, the shape of the initial 
configuration strongly affects the shape of the shock front. However, after some time 
the shock front sharpens to a discontinuity M irrespective of its initial shape. To 
understand the physical reason for this, consider yourself moving with the shock front 
and looking back towards the heated quark matter. Particles farther away from the 
front recede more slowly than particles just at the front. When the shock front is still 




(20) 
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smooth, it follows from energy- momentum conservation that entropy is conserved, i.e., 
with the quark matter equation of state, 

d t (T 3 1 )+d x (T 3 1 v)=0 . (21) 

This means that the temperature has to rise farther away from the front in order to 
accommodate the entropy of the matter moving with a lower velocity. However, as an- 
other consequence of energy-momentum conservation, temperature cannot rise enough 
to accommodate all the entropy inside a constant-sized volume, and a "traffic jam" 
-phenomenon occurs causing a discontinuity. An upper limit to the rate of jamming is 
clearly given by the difference of the velocity of the matter going into the jam and the 
velocity of the jam. This difference is just Vfl U id- Then the time scale of the sharpening 
is determined by ffl u id and on how smooth the shock front was in the beginning. The 
latter depends on how near the initial configuration was of the extremum bubble. Be- 
cause the shock front can initially be very wide and the flow velocity is very small, the 
time scale of the sharpening is very large. 

In fig. 10 a collision of two bubbles is shown. Due to the use of reflective 
boundary conditions, our grid corresponds to a situation in which several bubbles 
nucleate simultaneously at equal spacings, and therefore collisions are possible. The 
distance between the bubbles is in our picture Ax = 640 which is much less than the 
actual distances in the early universe, but this is not essential for the present analysis. 
At time t = 240, both the phase transition surface and the shock front are moving 
to the right. At t — 720, the shock fronts of neighbouring bubbles have collided, and 
the quark matter between the reflected shocks heats to a temperature higher than T c . 
At t = 960 the shock front and the phase transition surface are just about to collide; 
at t = 1040, the collision has just happened. The temperature of the hadron phase 
increases significantly. A shock front (a) continues to the hadron phase heating it up 
and making the matter move, but at the same time a rarefaction wave (b) is reflected 
back to the quark phase, cooling it down. The phase transition surface continues to 
move to the right, but its velocity has decreased from 0.06 to about 0.015. Some simple 
scenarios for the collision of a shock front and a deflagration front have been presented 
in ref. 0, and the course of events in fig. 10 is a combination of the scenarios a and c. 

As was seen in fig. 10, even one shock front can heat the hadronic matter consid- 
erably. This implies that only a few shock fronts are needed to raise the temperature 
of the hadronic matter to T c . This situation is shown in fig. 11, where several collisions 
are allowed to happen. The upper picture shows T(t,x) and the lower picture <fi(t,x). 
As the temperature of the hadron matter rises to T c the growth of the hadron bubbles 
is halted. In the cosmological context, a stage of slow growth of the hadron bubbles 
along with the expansion of the universe follows. In the electroweak case, the latent 
heat is smaller and the critical temperature is not reached. The bubbles fill the space 
and the fluid reaches the reheating temperature T re ^ eaut . However, we must remember 
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that the present simulations are 1+1 dimensional. In the 3-d case, the shock fronts are 
considerably weaker || and the time it takes for the universe to reheat is not necessar- 
ily the same as in the 1-d case. Neglecting the expansion of the universe, the fraction 
of space taken by the hadron matter immediately after reheating in the QCD case and 
the reheating temperature in the EW case are nevertheless the same as in the 1-d case. 

All the numerical solutions we have discussed so far were deflagrations. However, 
using a larger supercooling than eq. (|T8| ) indicates for the above (qh) values of a and 
L, we found detonation solutions. For Tf = 0.90 such a solution is shown in fig. 12. 
The upper picture shows T(t, x) and the lower picture 4>(t, x). The detonation expands 
with velocity t>det — 0.90 and a rarefaction front follows slightly behind. This is a weak 
detonation. 

As was seen in fig. 5, for some values of a and L even the physical nucleation 
temperature is so low that detonations could appear. In particular, even though the 
latent heat L is usually taken to be of the order of T c 4 , there is nothing in the lattice 
computations to rule out the possibility that L could as well be 0.1T C 4 or even less. Using 
the parameter values^ o = 0.1, L = 0.1, and £ c = 6, we did a sequence of computer 
runs varying the value of T (see fig. 13). This results in considerable supercooling, 
and the nucleation temperature (eq. (0)) is Tf = 0.891T C . For small values of V the 
bubbles grow as weak deflagrations. Increasing T increases t>d e f, until at about T = 10 
it approaches c s . Now t>d e f > c s would indicate a strong deflagration. Instead, as T is 
further increased, the solution shifts to a weak detonation. 



5 Steady-state variables of deflagration bubbles 

Very soon after the nucleation, the growing bubbles reach a steady-state configura- 
tion (see fig. 7). In the steady-state situation there exists a simple and very accurate 
way of finding out the interesting stationary variables T q , T h , t> s hock, ^fluid and t>d c f of 
the deflagration bubble apart from the above-presented integration of eqs. (|9|). This 
method can also be used to check the accuracy of the above integration. In the rest 
frame of the phase transition front the equations to be solved were given in eqs. ([TTHl3|) . 



Let us think that with the two conservation equations fll2|) and fllBf), we solve 
for T(x) and v (x) in terms of (f>(x) and 0'(x). Substituting these to eq. (|TT1), we get 
a second order differential equation for the field <p alone. However there are three 
boundary conditions: the derivative <fr'(x) must vanish asymptotically in both phases 
and we know the value of the field <f>(x) (only) in the quark phase. The system is over- 
determined, and only for certain values of the "constants" are there solutions. This is 
thus an eigenvalue problem. Assuming that we are given T g , we can solve for T h , t> def 



2 The larger value of £ c is needed for the potential of eq. (||) to be applicable for the range of 
temperatures in question. 
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and Wfluid- 

Next, consider the shock front. Because the shock front is in the quark phase, 
the field vanishes everywhere and we are left with very simple energy-momentum 
conservation equations. Solving them ||, §135] we get 



(22) 



Therefore, given Tf and T g , we can write down f s h oc k = Ui and t>fl ui d = (3/2) (^i — v 2 ). 

Now remember that a priori we only know Tf. However, guessing some T q , we 
get both from eqs. (p!THl3D and eq. (|22]) a value for Vfl U id- When we manage to guess such 



a T 9 that these two numbers agree the whole problem is solved. It is easy to make this 
method of solving the steady-state variables very accurate. Therefore we can use this 
method to check the accuracy of the dynamical integration with all time derivatives. 
With a reasonable number of grid points, the differences between the results of the two 
methods of integration on the steady-state variables are much less than 1% . Using 
this method, we can also easily calculate the steady-state variables as a function of T. 

5.1 Analytical approximations 

In this section we study what can be said analytically of the solutions of the steady- 
state equations (|TT|-|TB|). We use the bag equation of state ( ]IRO ) and assume that the 
velocities v g and Vh are non-relativistic (v 2 <C 1) and that the temperatures T q and Th 
are near T c (T q — Th <C T c ). The entropies are s q = 4:a q T g and Sh = 4a^T^ and the 



enthalpy is w = Ts. The one-parameter family of solutions of the two equations (12 



and (p~3| ) has been studied in detail in the literature [0, j|, ||, [12], [15], [L8|] , and the main 



problem is to find which solution the new equation (|TT|) picks out of this family. 

Evaluating eqs. (O) and ([13]) in the rest frame of the phase transition front for 
x = — oo and x = oo one obtains the usual energy- momentum conservation equations 

W q V q = W h V h 

(23) 

W q V l +Pq= W h vl + p h . 

From the upper equation it follows that 

vh = w 1= OqTl _ ^ r ^ (24) 
v q w h a h Tl a h 

Thus the fluid velocity is related to the deflagration front velocity by 

^fiuid = v h - v q = (1 - l/r)v dei . (25) 
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From the lower of eqs. ( p3|) we get 



Ph-Pg _ ahTh +L/A- a q T % 
w h - w g ^a h Tfi - Aa q T* 



= ^ = = — -. — =7 : — (26) 



'' " q Wh — w q 1,1 - r 1 

which in the limit T g , « T c = 1 implies that 

1 - T h = r(l - T 9 ) + r(r - l> g 2 . (27) 

To find the consequences of eq. (0), multiply it by (f>'(x) and integrate over the 
real axis. Using the equations <9U(0,T)/<90 = dVjd<\) - {dV/dT){dT/d<f)), V{-oo) = 
—L(l — Th), V(oo) = and 0'(±oo) = 0, and replacing the velocity v (x) by its absolute 
value (see fig. 8) one obtains 

L{1 -T h ) = - J_J<j>'{x)fv{x) dx - J °-^dT . (28) 

This is easy to solve in the limit r = Vh/v q — > 1 and Th,T q — > T c since then one can 
take v (x) ~ Vh ~ v q out of the first integral (which then gives the interface tension a) 
and neglect the last term. The result is 

v h nv q &v b = T-(l-T q ). (29) 
a 

This is the formula for bubble wall velocity derived earlier [II, TJ, [T3|, [HJ for the case of 
small change in flow velocity, which may be appropriate for the EW phase transition. 

In the more general case r > 1 note first that, because <j)'(x) is approximately 
symmetric around x = 0, the first term in eq. (|28|) can be approximated by (cr/T) • (v q + 
Vh)/'}- From fig. 8 one sees that at x = 0, dV/dT is about L/2 (in the hadron phase, 
dV/dT w L, since V{— oo) ~ — L(l — T/J, see above). We therefore approximate the 
second term by —n(L/2)(Tf l — T q ), where k is of order unity. Using eqs. ( j2"4|) and (p7|), 
eq. fl28|) then becomes 

(1 - K -)Lr{r - iy q - ^Av q + L(l - T q )[r + ~(1 - r)] = . (30) 

From this one can solve t>def = Vh = rv q as a function of T 9 . If further k — 1, eq. (|30|) 
simplifies into the form 

ttt) u ** " n^ def + r(1 _ Tq) = ' (31) 

which for small T (small t>d e f) again gives the result v q = v b in eq. (p9|). 

It is worth noticing that from eq. ((3^) we get a lower bound for T q . Namely, we 
know that eq. (BUf) has a solution and this gives 



<x 2 (r + l) 2 , . 

q ~ 16L 2 r 2 r(r-l)(l-K/2)[r + (K/2)(l-r)] ' 1 ' 
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The numerical value of this formula is useful for large T. However, even for moderate 
T it is essential to notice the existence of a lower bound: the temperature T q is not a 
free physical parameter, and when the whole expanding physical bubble is considered, 
the shock wave always heats quark matter just enough to reach the safe T q area. 

To obtain the result for a true deflagration bubble one finally has to eliminate 
the temperature T q from eq. ( |3~1~D by using eq. fl22|). Expanding v&ai& in powers of 1 — T q 
and 1 — Tf one gets 

Ufluid = >/3[(l-7»-(l-T,)] . (33) 
Using eq. (p5|) and substituting 1 — T q in terms of v d e f into eq. fl3l|) gives 



..r+l) def \LT ' v/3 , 
From this the final equation for v<m as a function of T is 




Vdei + r(l - T 



. 



(34) 



fdef 



a" r — 1 

if + ~\ 




4r(l - T 



r + 1 



r + 1 



(35) 



This is an excellent approximation for non-relativistic deflagration front velocities as 
will be seen in the next section (see fig. 15). In the limit r — > eq. (^) correctly 
reduces to the further approximation t>d e f ~ rT(L/a)(l — Tf). 



5.2 Numerical results for steady-state walls 
5.2.1 The quark-hadron phase transition 

In fig. 14 the temperatures Tf, T q and T^ are shown as a function of T for qh parameters. 
Small T means large friction and small velocities; large T means small friction and large 
velocities. When the deflagration front velocity is small, equation ( pop tells us that the 
fluid velocity is very small (for our present parameters, r « 1.1 so that i>fl U id « O.l^def)- 
Then from eqs. (|22|) we learn that T q has to be very close to Tf. For large V the 
situation is opposite: the velocities are larger and T q is higher. The temperatures T q 
and T h satisfy the relation (p7|) very accurately. Notice that entropy production at 
the shock front requires the condition T q > Tf and obviously one also has to obey the 
condition Th < T c , but nothing prevents T q from exceeding T c . This fact has been 
noticed before (see e.g. ref. |7[]), but usually it is not taken seriously. One reason may 
be that some authors neglect proper boundary conditions and thus confuse T q with Tf. 
Another is that it has been argued in refs. that a transition front between quark 

matter at the temperature T q > T c and hadron matter at the temperature Th < T c is 
impossible, basically because such a front would be mechanically unstable. However, 
both arguments are based on the assumption that at T c there exists a homogeneous 
mixed phase, and that in the transition zone between the quark phase and the hadron 



14 



phase, the state of the matter is at some point just in this homogeneous mixed phase. 
Then the transition front would be equal to two transition fronts, one from the quark 
phase to the mixed phase and the other from the mixed phase to the hadron phase. If 
there is a microscopic order parameter field — like in our model — the order parameter 
interpolates between the two minima of the effective potential in the transition zone, 
and there is no homogeneous mixed phase. The only mechanism by which a phase 
transition surface of this kind could in principle split is that a rarefaction wave detaches 
from it, and the analysis in ref. does not apply to rarefaction waves. Hence, there 



seems to be no reason why T q could not exceed T c , if the phase transition effectively 
includes an order parameter field. 

In fig. 15 the propagation velocity of the phase transition surface v^ e i is shown 
with solid line. With the dashed line we have drawn the deflagration front velocity 



from eq. fl35l) . This equation is seen to hold very well when the velocity v^ef is non- 
relativistic. The simple small-velocity approximation v^f ~ rT(L/ a)(l — Tf) is drawn 
with dotted line. In fig. 16 the fluid velocity t> fluid is shown as a function of T (solid 
line) and in fig. 17 the shock velocity v shock is drawn (solid line). The shock velocity is 
compared to the sound velocity. 

In fig. 16, the dotted line shows the entropy production. By entropy production 
we mean the relative change of the total entropy of a fluid element as the shock front 
and the phase transition surface sweep over it. One must note that the volume of the 
fluid element changes in the course of the process by the relative amount (see fig. 18) 

y&nal 1 — ^fluidA'shock ,„„x 

Pv = tj = 1 7 > 1 • ( 36 ) 

^initial J- — "Wfluid/ Vdei 

Therefore we define the entropy production by 

As es !^%Z_£ig>) . (37) 

In the quark-hadron phase transition where there is a conserved baryon number, this 
is just the relative change of entropy per baryon. The entropy production in eq. Q3"7|) 



is related to the quantity As pw = ShJh v h ~ s qlq v q (measured in the rest frame of the 
phase transition surface), which is often (e.g. in ref. |12| ) used to describe entropy 
production, and to the analogously defined quantity As s h measured in the rest frame 
of the shock front, by the equation 

As = ( + ^^py) /s q (T f ) . (38) 

The first term in the numerator is vanishingly small in comparison to the second. 

According to ref. |H| the quantity r\ = — T c (dv& e i / dT q )v& e { determines the sta- 



bility of expanding bubbles. If 77 > 1, the bubbles are stable at all length scales; if 
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rj < 1, large scale fluctuations are unstable. For our qh parameters, the quantity r] is 
drawn in fig. 17 with dotted line. While the analysis of ref. [[HJ is not suited for large 



T where T q can exceed T c , we notice, however, that for V < 0.83 our numerical results 
imply 7] to be less than unity and therefore large scale fluctuations should be unstable 
in that case. It would be interesting to expand the present code to include more space 
dimensions to see whether the expanding bubbles remain stable. 

5.2.2 The electroweak phase transition 

In figs. 19-20 we address the same questions as above but for the ew parameters. 
Qualitatively the behaviour is rather similar to what it was for the qh case. However, 
some differences worth a comment exist. The most important difference is that the 
electroweak phase transition is essentially a massless transition: the number of effective 
degrees of freedom changes very little at T c . Quantitatively this is expressed by the fact 
that r = a/(a — L/4) = 1.0018 is much closer to unity than in the qh case. Because of 
the smallness of r, there is a new temperature relevant for the phase transition. This 
is the reheating temperature T reheat (in the abrupt reheating scenario), defined by the 
equation e q (Tf) = e/i(T re h ea t)- For our present ew parameters, T re heat = 0.9965. 

Consider the temperatures (fig. 19) and the deflagration front velocity (fig. 20). 
Because the energy liberated in the phase transition is travelling in the compressed 
region behind the shock front, the temperature can never reach T re h ea t- This keeps 
Th low even for large V. On the other hand, it follows from eq. (^) that for large V the 
velocities t>d e f an d Vfluid grow large (now Vfl U id ~ 0.002t>d e f). Then it is seen from eqs. ( |2"2"| ) 
that the temperature T q has to rise considerably when T is large. But with r close to 
unity and the difference between T q and T h large, it is seen from eq. (|27D that the 
deflagration front velocity has to be very large indeed for large V. In fact the velocity 



Vdef becomes moderately relativistic and eqs. (p23H35|) are no longer strictly applicable. 
Another way to understand the connection between the high temperature T q and the 
large velocity t>d e f is to notice that the closer v^i is to the shock front velocity, the 
thinner is the area between the shock front and the phase transition surface. Therefore 
this thin area has a high temperature in order to accommodate all the latent heat 
released. 

Comparing figs. 14 and 19 we notice that for the ew parameters the values of V 



where t>d e f changes rapidly are larger than for the qh case. From eq. (35) we see that 
the relevant scale for V is roughly cr(r + 1)/L(r — 1). For the qh case the numerical 
value of this quantity is 1.1 and for the ew case 96 which explains the difference. 

Finally, let us compare our results to those of some other authors. Recently, the 
velocity of growing deflagration bubbles in the EW phase transition has been estimated 



for instance in refs. [[14], [Uj. In ref. |16| the authors note that the relation t>d c f > 1/30 
has to be satisfied in order not to diffuse away the baryon asymmetry and that probably 
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7t> ~ 1 — 2, that is, t> ~ 0.7 — 0.9. This corresponds to strong deflagrations, which seem 
to be an unlikely mechanism for bubble growth 
0.1 
r ft 



||. In ref. [14 velocities of order 
- 0.2 are obtained. Using simple kinetic theory to estimate the value of T we get 
10 



100 

very natural. 



11, 34], and then from fig. 20 mildly relativistic velocities appear to be 



6 Conclusions 

We have presented a model for phase transition bubbles in the early universe. In our 
model an order parameter field <fi with an effective potential V(<p, T) is coupled to 
a fluid with a dissipative constant T. Starting from an initial condition of a newly 
nucleated bubble we have numerically evolved the coupled hydrodynamical and field 
equations to follow the growth of the bubble in 1+1 dimensions. After some time the 
bubble reaches a stationary (similarity) state, where it grows at a constant velocity. We 
have then also studied the solutions to the corresponding stationary equations, both 
numerically and in analytical approximations. 

Typically the bubbles grow as weak deflagrations, and therefore with a subsonic 
velocity. The growth velocity is determined by the value of T, a large T (a weak coupling 
between and the fluid) leading to a large velocity and vice versa. Thus, we have been 
able to reduce the calculation of the growth velocity to the microscopic calculation of T. 
Our results show that the preheating caused by the shock front plays an essential role 
in the growth process, and that the temperature T q could even exceed T c . Reheating 
caused by collisions of expanding bubbles was also explicitly computed. 

If one uses simple kinetic theory to estimate T for the EW transition, one is 
naturally led to mildly relativistic velocities. For the QCD transition, one would di- 
mensionally expect that r m 1 [I5|. Then from fig. 15 for L = 2T C 4 , a = 0.1T C 3 we 



note that t>d c f ~ 0.06. The expectation thus is that the velocities in the QCD case are 
smaller than in the EW case. 

In some regions of the parameter space the solutions switch from weak deflagra- 
tions to weak detonations, as T is increased. It has been speculated that instabilities 
could turn expanding deflagration bubbles into detonations |17j]. We have now found 



that in some cases the bubbles could expand as detonations even from the beginning. 
Often it has been assumed that detonations in these phase transitions would have to 
be Jouguet detonations || [7|, as in chemical burning ||, but this appears not to 
be the case ||46|| . In the cosmological context, weak detonations require less extreme 
conditions than Jouguet detonations, making detonations more likely. 

For the QCD phase transition, there is an interesting parameter region that 
cannot presently be ruled out. Here the latent heat is rather small, say L = 0.1T C 4 , 
and the surface tension is similar, say a = 0.1T C 3 . Then the average distance between 
the nucleation centers would be of the order of 10 meters and large scale hadronic 
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inhomogeneities would result. With our model we find that such bubbles could grow 
as detonations, leading to a picture of the phase transition that is rather different from 
the usual one. 
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Figure captions 



Fig. 1: The velocities in the rest frames of the phase transition and shock fronts. 
Fig. 2: Structure of a deflagration bubble. 
Fig. 3: Structure of a detonation bubble. 

Fig. 4: The values of incoming (q) and outgoing (h) energy densities at the phase 
transition front, corresponding to deflagrations (upper left triangular region) and det- 
onations (lower right). This figure is for a bag model with r = a q /a h = 2. Point A 
corresponds to T q = T h = T c . In deflagrations the fluid has been heated by the shock, 
so T q > Tf, whereas for detonations T q = Tf. For a given initial temperature T/, there 
is a one-dimensional space of solutions, indicated by the dashed lines. The require- 
ment of non-negative entropy production restricts the solution below the line AS* = 0. 
Thus, in this case, only deflagrations are allowed for Tf = 0.9T C . For Tf = 0.8T C , both 
deflagrations and weak detonations are possible. For Tf = 0.7T C , Jouguet detonations 
are allowed, too. 

Fig. 5: For sufficiently small latent heat L, or sufficiently large surface tension a, 
detonations are possible in addition to deflagrations. The precise regions in the (a, L) 
plane depend on the equation of state. This figure is for a bag model with a q = 
51.257r 2 /90, ah = a q — L/4T^ (solid boundaries), or with a q = ah + L/4T^, ah = 
17.257r 2 /90 (dashed boundaries). 

Fig. 6: The meaning of the reflective boundary conditions. 

Fig. 7: The early stages of bubble growth. The upper picture shows T(t, x) and the 
lower picture shows v (t, x) for the qh parameters. 

Fig. 8: The variables <f>(x), T(x), v(x) and dV/dT as measured in the rest frame of 
the stationary phase transition surface. Hadron phase is on the left and quark phase 
on the right. By the velocity v we denote here the absolute value of the velocity: the 
direction of flow is from the quark to the hadron phase. 

Fig. 9: The sharpening of the shock front. 

Fig. 10: The collision of two bubbles. See the body of text for explanation. 
Fig. 11: Final stages of the phase transition for the qh parameters. 
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Fig. 12: A detonation solution for the qh parameters with T = 2 and Tf = 0.90. 
Note that m i n depends on temperature. 

Fig. 13: The bubble growth velocity for a = 0.1, L = 0.1 and £ c = 6 as a function 
of the parameter T. At approximately T = 10 the solution changes from a weak 
deflagration to a weak detonation. 

Fig. 14: The temperatures Tf, T c , T g , as a function of T for the qh parameters. 

Fig. 15: The velocity v^cf as a function of T for the qh parameters. The dashed line 
is the approximation from eq. fl35|) and the dotted line is the further approximation 
v Ae{ ^rT{L/a){l-T f ). 

Fig. 16: The fluid velocity and the entropy production for the qh parameters. 

Fig. 17: The shock velocity and the parameter r\ for the qh parameters. 

Fig. 18: The dashed lines indicate the average fluid flow in a similarity deflagra- 
tion solution. The volume of a fluid element changes when the shock front and the 
deflagration front sweep over it. 

Fig. 19: The temperatures Tf, T c , T q , Th and T re h ea t for the ew parameters. 

Fig. 20: The deflagration front velocity, the fluid velocity and the entropy production 
for the ew parameters. 
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